This is the documentation of the \texttt{MODULE SpecialFunc}, a set
of \texttt{FORTRAN 90} routines to compute the value of some
functions. This module make use of the \texttt{MODULE NumTypes},
\texttt{MODULE Constants}, \texttt{MODULE Error} so please read the
documentation of these modules \emph{before} reading this.

\section{Function \texttt{GammaLn(X)}}
\index{GammaLn@Function \texttt{GammaLn(X)}}

\subsection{Description}

Compute $\log(\Gamma(X))$.

\subsection{Arguments}

\begin{description}
\item[\texttt{X}:] Double (DP) precision. The point in which we want to
  know the value of $\Gamma(X)$.
\end{description}

\subsection{Output}

A real Double precision (DP).

\subsection{Examples}

This program should write the factorial of the first $100$ numbers.

\begin{lstlisting}[emph=GammaLn,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Gamma Function.,
                   label=gammaln]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Integer :: q


  Do q = 1, 100
     Write(*,'(1A13,1I4,1A3,1ES33.25)')'Factorial of:', q, ' = ',&
          & exp(GammaLn(Real(q+1,kind=DP)))
  End Do


  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{erf(X)}}
\index{erf@Function \texttt{erf(X)}}

\subsection{Description}

Compute values of the error function $\text{erf}(X)$. A high
precision algorithm by W. J. Candy~\cite{cody:1969aa} is used.

\subsection{Arguments}

\begin{description}
\item[\texttt{X}:] Single (SP) or double (DP) precision. The point in
  which we want to know the value of $\text{erf}(X)$.
\end{description}

\subsection{Output}

Same type as input.

\subsection{Examples}

This program should compute the $\text{erf}(X)$ and its complementary
in an given point.

\begin{lstlisting}[emph=erf,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Error Function.,
                   label=erf]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Real (kind=DP) :: X

  Write(*,*)'Enter a number: '
  Read(*,*)X
  Write(*,'(ES33.25)')erf(X)
  Write(*,'(ES33.25)')erfc(X)

  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{erfc(X)}}
\index{erfc@Function \texttt{erfc(X)}}

\subsection{Description}

Compute values of the complementary error function $\text{erf}(X)$. A high
precision algorithm by W. J. Candy~\cite{cody:1969aa} is used.

\subsection{Arguments}

\begin{description}
\item[\texttt{X}:] Single (SP) or double (DP) precision. The point in
  which we want to know the value of $\text{erfc}(X)$.
\end{description}

\subsection{Output}

Same type as input.

\subsection{Examples}

This program should compute the $\text{erf}(X)$ and its complementary
in an given point.

\begin{lstlisting}[emph=erfc,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the complementary Error Function.,
                   label=erfc]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Real (kind=DP) :: X

  Write(*,*)'Enter a number: '
  Read(*,*)X
  Write(*,'(ES33.25)')erf(X)
  Write(*,'(ES33.25)')erfc(X)

  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{inverf(X)}}
\index{inverf@Function \texttt{inverf(X)}}

\subsection{Description}

Compute values of the inverse error function $\text{inverf}(X)$. A
combination of a rational approximation by Peter John Acklam, and
Halley's rational method is used to attain double precision.

\subsection{Arguments}

\begin{description}
\item[\texttt{X}:] Single (SP) or double (DP) precision. The point in
  which we want to know the value of $\text{inverf}(X)$.
\end{description}

\subsection{Output}

Same type as input.

\subsection{Examples}

This program should use the $\text{inverf}(X)$ function and test it is
ok.

\begin{lstlisting}[emph=inverf,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Inverse Error Function.,
                   label=inverf]
Program Test

  USE NumTypes
  USE SpecialFunc

  Real (kind=DP) :: X = 0.6783457843_DP, Y

  Write(*,'(ES33.25)')inverf(-0.998_DP)
  Write(*,'(ES33.25)')inverf(-0.5_DP)
  Write(*,'(ES33.25)')inverf(0.988_DP)

  Y = inverf(X)
  Write(*,*)'Same number two times!'
  Write(*,'(ES33.25)')X
  Write(*,'(ES33.25)')erf(Y)

  Stop
End Program Test
\end{lstlisting}

\section{Function \texttt{Theta(i, z, tau[, Prec])}}
\index{Theta@Function \texttt{Theta(i, z, tau[, Prec])}}

\subsection{Description}

Compute the value of the $i^{\underline{\text{th}}}$ Jacobi theta
function $(i=1,2,3,4)$ with nome $q=e^{i\pi\tau}$ 
\begin{equation}
  \vartheta_i(z|\tau)
\end{equation}
For a definition and properties of these functions take a
look~\cite{ww:analysis}, here we will only say that following the
conventions of the cited reference, our Theta functions have
quasi-periods $\pi$ and $\tau\pi$. 

\subsection{Arguments}

\begin{description}
\item[\texttt{i}:] Integer. Which theta function we
  want to compute. $i$ must have one of the following values: $1,2,3,4$.
\item[\texttt{z}:] Complex Double Precision (DPC) or Complex Single
  Precision (SPC). The point in which we want to compute the Theta
  function.
\item[\texttt{tau}:] Complex, with the same precision as
  \texttt{z}. is the quasi period of the Theta function. must be in
  the upper half plane $(\textrm Im(\tau)> 0)$.
\item[\texttt{Prec}:] Real, Optional. If \texttt{z} is DPC (SPC),
  \texttt{Prec} must be double precision (single precision). An
  estimation of the desired precision of the result. The default value
  is $1\times 10^{-3}$
\end{description}

\subsection{Output}

If \texttt{z} is Double Precision Complex (SPC), the the result will be
Double Precision Complex (SPC).

\subsection{Examples}

\begin{lstlisting}[emph=Theta,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Jacobi Theta functions.,
                   label=theta]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Complex (DPC) :: Z, tau


  Z = Cmplx(0.546734, 2.76457643, kind=DPC)
  tau = Cmplx(0.0_DP, 3.76387540_DP)

  ! Check the quasi-periodicity of the Third
  ! Jacobi Theta function.
  Write(*,*)Theta(3, Z           , tau)
  Write(*,*)Theta(3, Z+Cmplx(PI_DP), tau)
  Write(*,*)Theta(3, Z+PI_DP*tau, tau) * &
       &exp(PI_IMAG_DPC*tau + 2.0_DP*UNITIMAG_DPC*Z)


  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{ThetaChar(a, b, z, tau[, Prec])}}
\index{Theta@Function \texttt{ThetaChar(a, b, z, tau[, Prec])}}

\subsection{Description}

Computes the value of the Theta function with Characteristics $(a,b)$
and quasi-periods $(\pi,\pi\tau)$ in the point $z$:
\begin{equation}
  \thetachar{a}{b}\left(z|\tau\right)
\end{equation}

\subsection{Arguments}

\begin{description}
\item[\texttt{a, b}:] Complex or Real, Single or double precision. The
  two characteristics of the Theta function.
\item[\texttt{z}:] Complex (Single or Double precision). The point in
  the complex plane.
\item[\texttt{tau}:] Complex (Single or Double precision). The
  quasi-period of the theta function. Must have $(\textrm{Im}(\tau)>0)$.
\item[\texttt{Prec}:] Real (Single or Double precision). Optional. An
  estimation of the desired precision of the value of the theta
  function. 
\end{description}

\subsection{Output}

Complex Single or Double precision, the same as the input values.

\subsection{Examples}

\begin{lstlisting}[emph=ThetaChar,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the Jacobi Theta functions with characteristics.,
                   label=thetachar]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Real(kind=DP) :: Deriv, X1, X2
  Complex (DPC) :: Wmas, Wmenos, Z, tau
  Integer :: q, s


  Z = Cmplx(0.546734, 2.76457643, kind=DPC)
  tau = Cmplx(0.0_DP, 3.76387540_DP)


  Write(*,*)'Theta 1:'
  Write(*,*)Theta(1, Z, tau)
  Write(*,*)-ThetaChar(0.5_DP,0.5_DP, Z, tau)
  Write(*,*)'Theta 2:'
  Write(*,*)Theta(2, Z, tau)
  Write(*,*)ThetaChar(0.5_DP,0.0_DP, Z, tau)
  Write(*,*)'Theta 3:'
  Write(*,*)Theta(3, Z, tau)
  Write(*,*)ThetaChar(0.0_DP,0.0_DP, Z, tau)
  Write(*,*)'Theta 4:'
  Write(*,*)Theta(4, Z, tau)
  Write(*,*)ThetaChar(0.0_DP,0.5_DP, Z, tau)
  

  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{Hermite(n,x[, Dval])}}
\index{Hermite@Function \texttt{Hermite(n,x[, Dval])}}

\subsection{Description}

Returns the value of the $n^{\underline{\text{th}}}$ Hermite
polynomial in the point $X$. If \texttt{Dval} is specified, the value
of the Derivative of the $n^{\underline{\text{th}}}$ Hermite
polynomial in the point $X$ is also returned.

\subsection{Arguments}

\begin{description}
\item[\texttt{n}:] Integer. Which Hermite polynomial wants to compute.
\item[\texttt{X}:] Real (Single or Double precision). The point in
  which we want to compute the Polynomial.
\item[\texttt{Dval}:] Real (Single or Double precision). Optional. If
  specified, it stores the value of the Derivative of the
  Polynomials.
\end{description}

\subsection{Output}

Real single or double precision (same as input). The value of the
$n^{\underline{\text{th}}}$ Hermite Polynomial in the point X.

\subsection{Examples}

\begin{lstlisting}[emph=Hermite,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the first 31 Hermite numbers.,
                   label=hermite]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Integer :: q


  Write(*,*)'The first 31 Hermite Numbers'
  Write(*,*)'http://www.research.att.com/~njas/sequences/A067994'
  Do q = 1, 31
     Write(*,'(1I4,1ES33.25)')q, Hermite(q, 0.0_DP)
  End Do

  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{HermiteFunc(n, x[, Dval])}}
\index{Hermite@Function \texttt{HermiteFunc(n, x[, Dval])}}

\subsection{Description}

Returns the value of the $n^{\underline{\text{th}}}$ Hermite
function
\begin{equation}
  \frac{1}{\sqrt{n!2^n\sqrt{\pi}}}e^{-x^2/2}H_n(x)
\end{equation}
in the point $X$. If \texttt{Dval} is specified, the value
of the Derivative of the $n^{\underline{\text{th}}}$ Hermite
function in the point $X$ is also returned.

\subsection{Arguments}

\begin{description}
\item[\texttt{n}:] Integer. Which Hermite function wants to compute.
\item[\texttt{X}:] Real (Single or Double precision). The point in
  which we want to compute the Polynomial.
\item[\texttt{Dval}:] Real (Single or Double precision). Optional. If
  specified, it stores the value of the Derivative of the
  function.
\end{description}

\subsection{Output}

Real single or double precision (same as input). The value of the
$n^{\underline{\text{th}}}$ Hermite function in the point X.

\subsection{Examples}

\begin{lstlisting}[emph=HermiteFunc,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Compute the Hermite functions.,
                   label=hermitefunc]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Real(kind=DP) :: Deriv, X1, X2, Sum
  Complex (DPC) :: Wmas, Wmenos, Z, tau
  Integer :: q, s


  Write(*,*)'A (really bad) proof of orthonormality:'
  X1 = -10.0_DP
  Sum = 0.0_DP
  Do q = -1000, 1000
     Sum = Sum + HermiteFunc(6,X1)**2
     X1 = X1 + 0.01_DP
  End Do

  Write(*,'(1ES33.25)')Sum*0.01_DP

  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{Basis(X1, X2, n, s, q, itau[, Prec]) }}
\index{Basis@Function \texttt{Basis(X1, X2, n, s, q, itau[, Prec]) }}

\subsection{Description}

Return the value of the basis elements of the Hilbert space $\mathcal
H_q$ of quasi-periodic functions
\begin{equation}
  \ket{n,s} = e^{\imath
    \frac{f}{2}x_1x_2} \sum_{k\in s + q\mathbb Z} 
  e^{-u^2/2}H_n(u)
  e^{2\pi\imath k\frac{x_1}{l_1}}\qquad n=0,\dots\infty;s=1,\dots,q
\end{equation}
defined in the appendix of~\cite{Gonzalez-Arroyo:2004xu} (look there
for more details and properties).

\subsection{Arguments}

\begin{description}
\item[\texttt{X1,X2}:] Real (Single or Double precision). The point in
  the Torus.
\item[\texttt{n,s}:] Integer. Specify which element of the basis.
\item[\texttt{q}:] Integer. Specify the Hilbert space $\mathcal H_q$.
\item[\texttt{itau}:] Real (Single or Double precision). Specify the
  ratio of quasi-periods: $\mathtt{itau} = l_2/l_1$.
\item[\texttt{Prec}:] Real (Single or Double precision). Optional. An
  estimation of the desired precision. 
\end{description}

\subsection{Output}

Complex single or double precision, depends of the input  arguments.

\subsection{Examples}

\begin{lstlisting}[emph=Basis,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the bassi of a special Hilbert
                   space (details in~\cite{Gonzalez-Arroyo:2004xu}).,
                   label=basis]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Real(kind=DP) :: X1, X2
  Complex (DPC) :: Wmas, Wmenos,
  Integer :: I, q, s


  Write(*,*)'Looking at the quasi-periodicity properties:'
  X1 = 0.97834D0
  X2 = 0.873873D0
  q = 4
  s = 3
  Do I = 0, 8
     Wmas   = Basis( X1,  X2+1.0_DP, I, s, q, 1.0_DP, 1.0D-15) * &
          & exp(PI_IMAG_DPC*X1*q)
     Wmenos = Basis( X1+1.0_DP,  X2, I, s, q, 1.0_DP, 1.0D-15) * &
          & exp(-PI_IMAG_DPC*X2*q)
     Write(*,'(1I3,2ES33.25)')I, Basis( X1,  X2, I, s, q, 1.0_DP, 1.0D-15)
     Write(*,'(1I3,2ES33.25)')I, Wmas
     Write(*,'(1I3,2ES33.25)')I, Wmenos
  End Do


  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{Factorial(N)}}
\index{Factorial@Function \texttt{Factorial(N)}}

\subsection{Description}

Compute $N!$. Better (faster and more accurate for small numbers) than
the use of \texttt{GammaLn} to compute the factorial of a number.

\subsection{Arguments}

\begin{description}
\item[\texttt{N}:] Integer. The number to compute the factorial.
\end{description}

\subsection{Output}

A real Double precision (DP).

\subsection{Examples}

This program should write the factorial of the first $100$ numbers.

\begin{lstlisting}[emph=Factorial,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing the factorial.,
                   label=gammaln]
Program TestSpecialFunc

  USE NumTypes
  USE SpecialFunc

  Integer :: q


  Do q = 1, 100
     Write(*,'(1A13,1I4,1A3,1ES33.25)')'Factorial of:', q, ' = ',&
          & Factorial(q)
  End Do


  Stop
End Program TestSpecialFunc
\end{lstlisting}

\section{Function \texttt{Legendre(l, m, X)}}
\index{Legendre@Function \texttt{Legendre(l, m, X)}}

\subsection{Description}

Computes the value at $X$ of the $(l,m)$ Legendre polynomial $P_l^m(X)$.

\subsection{Arguments}

\begin{description}
\item[\texttt{l,m}:] Integer. Specifies which Legendre polynomial we
  want to compute.
\item[\texttt{X}:] Real sinple or double precision. The point at which
  we want to compute the polinomial.
\end{description}

\subsection{Output}

A real simple or double precision (same as \texttt{X}).

\subsection{Examples}

This program writes the values of the Legendre polynomials
$P_l^m(0.5)$ for $0\le m \le l$.

\begin{lstlisting}[emph=Legendre,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing some legendre polynomials.,
                   label=legendre]
Program TestLeg


  USE NumTypes
  USE SpecialFunc

  Real (kind=DP) :: Y

  Read(*,*)L

  Do I = 0, L
     Y = Legendre(L, I, 0.50_DP)
     Write(*,*)I, Y, Legendre(L, I, 0.5_SP)
  End Do

  Stop
End Program TestLeg
\end{lstlisting}

\section{Function \texttt{SphericalHarmonics(l, m, th,ph)}}
\index{Spherical@Function \texttt{SphericalHarmonics(l, m, th, ph)}}

\subsection{Description}

Computes the value at $(\theta,\phi)$ of the $(l,m)$ spherical
harmonic $Y_l^m(\theta,\phi)$. 

\subsection{Arguments}

\begin{description}
\item[\texttt{l,m}:] Integer. Specifies which spherical harmonic we
  want to compute.
\item[\texttt{th,ph}:] Real single or double
  precision. $(\theta,\phi)$ coordinates of the point at which we
  we want to compute the vlue of the spherical harmonic.
\end{description}

\subsection{Output}

A complex simple or double precision (same as \texttt{th}).

\subsection{Examples}

This program checks the relation $Y_l^m(\theta,phi) =
(-)^mY_l^{*-m}(\theta,\phi)$ for $\theta=0.5, \phi=0.35$.

\begin{lstlisting}[emph=SphericalHarmonic,
                   emphstyle=\color{blue},
                   frame=trBL,
                   caption=Computing some spherical harmonics.,
                   label=sph]
Program TestSph


  USE NumTypes
  USE SpecialFunc

  Complex (kind=DPC) :: Y

  Read(*,*)L

  Do I = 0, L
     Y = SphericalHarmonic(L, I, 0.50_DP, 0.35_DP)
     Write(*,'(1A,1I4,1A,1I4)')'Spherical Harmonic l=', L, ' m=+-', I
     Write(*,*)'  +m ', Y
     Write(*,*)'  -m ', SphericalHarmonic(L, -I, 0.5_DP, 0.35_DP)
  End Do

  Stop
End Program TestSph
\end{lstlisting}


% Local Variables: 
% mode: latex
% TeX-master: "lib"
% End: 

